Non-universal level statistics in a chaotic quantum spin chain 
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We study the level statistics of an interacting multi-qubit system, namely the kicked Ising spin 
chain, in the regime of quantum chaos. Long range quasi-energy level statistics show effects analo- 
gous to the ones observed in semi-classical systems due to the presence of classical periodic orbits, 
while short range level statistics display perfect statistical agreement with random matrix theory. 
Even though our system possesses no classical limit, our result suggest existence of an important 
non-universal system specific behavior at short time scale, which clearly goes beyond finite size 
effects in random matrix theory. 
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I. INTRODUCTION 

One of the key discoveries of quantum chaos has been 
the so-called quantum chaos conjecture, originally pro- 
posed in Refs. 1}. It claims that even simple non- 
integrable quantum systems, whose dynamics is suffi- 
ciently complex (say, dynamically mixing in the classical 
limit), possess quantum fluctuations which can be de- 
scribed by a universal ensemble of random matrices with- 
out any free parameters [2j ■ Although a strict mathemat- 
ical proof of this conjecture is still missing, its theoretical 
understanding has recently been considerably deepened 
0. Still, it is known from the early years of quantum 
chaos [H, that level fluctuations exhibit universal fea- 
tures only on sufficiently small energy scales, or long time 
scales, whereas one obtains system specific non- universal 
features on long energy scales (short time-scales) which 
can be usually understood and computed in terms of clas- 
sical orbits. 

Therefore it seems that the picture is quite complete 
and satisfactory for systems possessing a well defined 
classical limit. But what about simple systems which 
do not have a classical limit, e.g., systems of interacting 
fcrmions, or systems of interacting qubits? In such sys- 
tems, dynamical complexity can be reached in the ther- 
modynamic limit of many interacting particles In 
some exactly solvable cases formal similarities between 
the thermodynamic limit and the semi-classical limit can 
be established [f|. For example, one may start by consid- 
ering simple, non-integrable, many-particle Hamiltonians 
with local interaction which are specified by only a few 
(non-random) parameters. Can quantum spectral fluc- 
tuations of such systems be described by universal en- 
sembles of random matrices? If yes, what are the energy 
scales of such universality? Is there a breaking of univer- 



sality at sufficiently large energy ranges? How does the 
universality breaking scale in the thermodynamic limit? 
In this paper we address these questions in a simple dy- 
namical system, namely an Ising chain of spin 1/2 parti- 
cles on a Id ring, kicked periodically with a homogeneous, 
tilted magnetic field. We performed careful numerical 
calculations of quasi-energy spectra and their statisti- 
cal analysis. For appropriate values of model's param- 
eters, corresponding to strong integrability breaking, we 
indeed find both, the universality regime for sufficiently 
small energy scales, where no statistically significant de- 
viations from random matrix prediction of infinitely di- 
mensional Circular Orthogonal Ensemble (COE) could 
be detected, and a non-universality regime for large en- 
ergy scales (or small times, corresponding to few kicks), 
where clear, statistically significant deviations from Ran- 
dom Matrix Theory (RMT) prediction have been found. 
Most notably, our analysis shows that the spectral form 
factor exhibits significant deviations from RMT at the 
time scale corresponding to one or few kicks (Floquet pe- 
riods). This result could be intuitively understood as a 
qubit analogy of "shortest periodic orbit" correction, but 
its precise theoretical understanding is at present open. 



II. THE SYSTEM 

The system we study is a kicked Ising chain (KIC) 0] , 
namely a ring of L spin 1/2 particles which interact with 
their nearest neighbors via a homogeneous Ising interac- 
tion of dimensionless strength J and being periodically 
kicked with a homogeneous magnetic field of dimension- 
less strength b. During the free evolution, i.e. between 
the kicks, the system evolves with the unitary propagator 
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and the action of the kick is described by the unitary 
operator 

I ^ \ 

£4ick0) = exp \-i 22 b ■ a j J , (2) 

with a x ^ v,z being the Pauli matrices of particle j and 
(ij = (<jj,(jj,(jj). The Floquet operator for one period 
is thus 

Uki = U lsins (J)U kick (b). (3) 

We must also impose periodic conditions in order to close 
the ring: <tl = &o- During this paper we shall use the 
so called computational basis, which is composed of joint 
eigenstates of <j| . This set of basis states can be written 
as S = {\momi . . . rriL-i), withm^ G {0, 1}}. 

In order to understand the spectrum one must first 
discuss the symmetries in the system. Let us start 
with the translational symmetry. The corresponding 
operator T is defined on the computational basis as 
T|momi • • • rriL-i) = |rn£_imo • • • m£_2), and is ex- 
tended to the entire Hilbert space by linearity. The ac- 
tion of this operator is to rotate the particles in the ring 
by one site. The eigenvalues of T are exp(27nfc/L) with 
k G Ii/l- Hence the Hilbert space is foliated into L sub- 
spaces Ti — (J) fcgZ/i Ti-k- The evaluation of the dimen- 
sionality of each of these subspaces is described in Ap- 
pendix [X] The evolution operator © is translationally 
invariant, and hence [Uki, T] — 0. 

The next symmetry is an external reflection R. Its 
action on the basis S reads R\momi ■ ■ ■ = 
|m£^i?Tt£_2 • ■ • mo), and also [Uki, R] = 0. The two sym- 
metries, T and R, do not commute. It must be noticed 
that if G Hk, then R\tp) G H-k- This provides an 
additional symmetry within the subspace Hq (and H.l/2 
for even L) . Thus these marginal subspaces are regarded 
as "special" and have slightly different properties than 
the rest. We shall not consider them for the purpose of 
statistical analysis in this article. 

Finally, we define the anti-unitary symmetry K! '. It 
acts as a mirror reflection within each spin with respect 
to the plane that contains both b and the unitary vector 
in the z direction (the direction of the Ising interaction) . 
We can rotate our coordinate system around the z axis 
in each qubit so that b only has components in the x 
and z directions. Then, K is simply complex conjuga- 
tion, provided that a x and a z are set real, as is the usual 
choice. This symmetry operation also changes the sign 
of the momentum. Composing /C = IC'R, we arrive to 
an anti-unitary symmetry that preserves the momentum 
i.e. an anti-unitary symmetry within each TLk- 

Concluding, Uki has a rotational symmetry R that 
foliates the spectrum into L different sectors; the sec- 
tor k has identical spectrum as the sector —k. Hence 
for a fixed number of qubits we expect to have a maxi- 
mum of (L — l)/2 relevant sectors, each sector having a 
Hilbert space dimension M » 2 L /L. Since each sector 
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FIG. 1: (Color online) We study the behavior of the in- 
tegrated nearest neighbor spacing distribution, W(s). The 
noisy curve shows the difference between the numerical data 
for 18 qubits, averaged over the different relevant Ttk spaces, 
and the Wigner surmise. The smooth (red) curve is the dif- 
ference between infinitely dimensional COE solution and the 
Wigner surmise. The expected standard deviation due to the 
finite size of the spectrum aw [see Eq. Q] is also indicated 
as the shaded area surrounding the RMT result. In the inset 
we present a similar figure with the results for each of the TLk 
subspaces plotted separately, together with the error associ- 
ated with each individual spectrum, as the thick black curves. 
(Different colors represent different sectors, according to the 
coding shown in Fig. [3]) 



has an internal anti-unitary symmetry, we shall compare 
the statistical properties of system's spectrum to those 
of Dyson's Circular Orthogonal Ensemble [2] of random 
matrices, of appropriate dimension. 

For the rest of the presentation we fix parameter val- 
ues of our system J = 0.7, b — (0.9, 0, 0.9) for which the 
integrability of the model is strongly broken. We believe 
that for these parameter values the system is a generic 
representative of quantum chaos. Using highly optimized 
numerical methods (see appendix [Bj we have been able 
to diagonalize the model accurately for sizes up to L = 18 
qubits. The eigenvalues of the Floquet propagator Uki 
have been written as exp(— up n ), where Lp n are known 
as quasi-energies, and have been grouped with respect 
to the known quasi-momentum k. Statistical analyses of 
desymmetrized quasi-energy spectra {^p n } and their in- 
terpretation are given in the following sections. In order 
to compare with the RMT formulae we normalize the 
quasi-energies, i.e. write s n = in order to have 

mean level spacing equal to one (AT denotes the dimen- 
sion of the Hilbert space). Similarly, s = will denote 
the spectral variable, unfolded to a unit mean level spac- 
ing. 

III. UNIVERSALITY REGIME 

Let us first analyze the most commonly studied spec- 
tral statistics of chaotic systems, that is the nearest 
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FIG. 2: (Color online) In this plot we show the behavior of 
the form factor for 18 qubits. The black dots show its value at 
integer times for the quasi- momentum sector k = 1. In order 
to appreciate clearly its behavior it is necessary to perform 
a windowing over short ranges of time (th/25). The results 
for each of the fc-spaces are shown according to the symbol 
scheme in Fig. [3] The average over the different spaces as 
well as the theoretical curve is plotted as a black and red line, 
respectively. In order to compare with the ensemble fluctua- 
tions, we plot in the inset the difference from the theoretical 
prediction of both, the spectra for the KIC (in black) , and the 
spectra of an equal number of random realizations of COE 
members with the same dimension (in red). 
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FIG. 3: (Color online) We observe the variance E 2 for each 
symmetry sector (as different symbols) and its average (thick 
black curve), for 18 qubits. Note that lower abscissa indicates 
the unfolded spectral variable, while the upper abscissa the 
one without unfolding. The average curve is almost indistin- 
guishable from the theoretical value (thin red curve). In the 
inset we compare the deviation of the averag ed E 2 for both 
the KIC (in black) and the COE (in red) from the theoretical 
value Eqoe- No qualitative difference is observed. 



and for infinitely dimensional COE has the form 



neighbor level spacing distribution P(s). P(s)ds is the 
probability that the distance between two nearby (un- 
folded) quasi- energies is between s and s + ds. P{s) 
has been computed for the KIC and compared to the 
exact random matrix COE result (computed from Pade 
approximants |8j) with satisfactory results (not shown). 
However, since the details of such plot depend on the size 
of the binning of histograms, we prefer to show the cu- 
mulative (integrated) level spacing distribution W(s) = 
Jq ds'P(s'). In fig. Q] we show a comparison of W(s), 
both for the KIC and the exact infinitely-dimensional 
COE result, with the Wigner surmise Wwigncr(s) = 
1 — exp(— 7rs 2 /4). The expected statistical fluctuation 
of cumulative probability can be estimated 0] as 



aw 



W(l - W) 



(4) 



and gives a very realistic estimate of actual fluctuations 
of our dynamical system. We plot the results both for 
individual quasi-momentum k subspaces, and averaged 
over all k. In conclusion, based on the nearest neighbor 
level spacing distribution we find no significant deviations 
from universality, i.e. from COE model statistics. 

Further on, we have studied other statistical measures 
of quasi-energy spectra, which are more sensitive to long- 
range spectral correlations, namely the number variance 
and the form factor 0. The spectral form factor K 2 , is 
defined for discrete time t as K 2 (t/r H ) = |Tr [/*| 2 /-A/"> 
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th = TV denotes the discrete Heisenberg time, namely 
the number of kicks in which the average quasi-energy 
level separation grows to 27r. Asymptotic finite dimen- 
sion corrections to the form factor have been computed, 
and for small r <C 1, the result reads 
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The number variance £ 2 (s) gives the variance of the num- 
ber of levels in an unfolded spectral interval of length 
s. The RMT formula for infinitely dimensional COE 
predicts a monotonically increasing variance Sq OE (s) = 
(2/7T 2 )[ln(27rs) + 1 + 7 - ^ 2 /8] + CO" 1 ), where 7 = 
0.5772 ... is the Euler constant Q. However, for a finite 
spectrum of Af quasi-energy levels this is not possible, 
since when the energy difference reaches the range of the 
spectrum the number of levels counted will always be the 
maximum and hence the number variance will be zero. 
For arbit rary finite dimension A/", there is an exact rela- 
tionship [n| between E 2 and K 2 that accounts for the 
finite range of the spectrum: 
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Truncating the above series at finite m with the form 
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FIG. 4: We plot E 2 for the KIC for L — 14 (upper plot) and 
L = 18 (lower plot). Different sectors are plotted using thin 
colored curves (with the same color coding as in Fig. [3j and 
the average value as the thick black curve. The theoretical 
(COE) prediction for Af dimensional circular random matrices 
[eq. 0] is also plotted as a smooth red curve. We observe 
saturation of the stiffness for the KIC, characteristic of semi- 
classical systems. 

factor given by eq. ^ provides an excellent asymptotic 
approximation to the COE number variance for finite Af. 

In figure [2] we compare the spectral form factor of the 
KIC with the infinitely dimensional COE, on a global 
time scale (on the order of Heisenberg time th = Af) . Of 
course, since the form factor is not self-averaging we had 
to perform some averaging over short time-windows in 
order to wash away the statistical fluctuations. We find 
no notable deviation from the COE. In order to estimate 
the expected fluctuations due to a finite sample of sys- 
tems (namely a set of ~ L/2 quasi- momenta k) we have 
also generated a similar average over the same number 
of random matrices of equivalent size. In the inset we 
plot the deviations of the form factor computed for the 
KIC, and the corresponding finite average over random 
members of the COE, from the exact RMT prediction. 
We observe that both behave similarly. In addition, we 
find very good agreement for the number variance E 2 of 
the KIC with the infinitely dimensional COE on short 
and intermediate spectral ranges s < 10 (see fig. [3]). In 
the inset the finite size fluctuations are also compared 
with the ones computed from appropriate finite samples 
of finite dimensional COE. Again we observe agreement. 

IV. DEVIATIONS FROM UNIVERSALITY 

As explained in the previous section, we expect that 
the number variance for a finite spectrum reaches a max- 
imal value at <p = tt (i.e. s ~ Af/2). Actually as already 
mentioned, one can compute a good analytical approx- 
imation to COE averages of number variance for finite 
Af using eq. (|7|), and the saturation can be understood 
as a consequence of discreteness of time in the sine-like 
transformation on the RHS of Q. 
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FIG. 5: (Color online) We see the dependence of E 2 *, as a 
function of the number of qubits L, both for the KIC and the 
COE of appropriate dimension. In the inset we plot ifioo. The 
constant line </? c = 0.39 is shown in red, as well the plot of 
£coe(A/V>c/27t, AT) in the main panel. 



Do the spectra of KIC in the regime of quantum chaos 
follow the same saturation as would be expected for typ- 
ical members of the COE or not? We have performed de- 
tailed numerical checks of these questions and report the 
results in the following figures. In fig. [4] we plot the num- 
ber variance for two different number of qubits (14 and 
18), for the KIC. We find a very clear and notable differ- 
ence: the data for the KIC tends to saturate at different, 
lower value of the unfolded spectral parameter <C Af, 
than COE, which typically saturate only at s ~ Af/2. 
Furthermore, the plateau is quite notorious. In the next 
plot (fig. [5]) we have determined the saturation threshold 
Soo and the saturation value E 2 ^ = E 2 (soo) as a function 
of the number of qubits L. 

Numerical results suggest that ss 0.062A/" oc 2 L /L, 
namely that Sqo is proportional to Af though it is smaller 
by a large constant factor. The saturation value E 2 ^ thus 
increases logarithmically with Af, or linearly with L. 

Perhaps a more clear picture is obtained after going 
into the time domain and inspecting the form factor 
Ki for a few kicks, which correspond to large spectral 
ranges of E 2 . This regime is analogous to the non- 
universality regime corresponding to the shortest clas- 
sical periodic orbit in quantum chaotic systems with 
well defined (semi) classical limit. However, we should 
not forget that our spin chain does not have any well 
defined classical limit or semi-classical regime. Still, 
it seems that ^(I/th), K 2 (2/th), etc., notably de- 
viate from expectations of COE of the same dimen- 
sions as the KIC propagator. Indeed, in figEl we show 
^(I/th) as a function of L, both for individual quasi- 
momentum k subspaces and the average over all rel- 
evant k, and find very clear and systematic deviation 
from COE expectation -K2,coe(1/ t h) = 2/Af. Further- 
more, deviation of ^(I/th) is much bigger that expected 
COE fluctuation of ^(Ith) which can be computed as 
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FIG. 6: We plot the value of K% evaluated at the first kick, 
for different quasi-momentum sectors (symbols according to 
fig. 0) and averaged over all quasi-momenta (filled circles). 
The red line indicates the theoretical COE value surrounded 
by one expected standard fluctuation (according to theoreti- 
cal COE fluctuation) indicated by gray area. The average is 
systematically bellow the expected RMT value. 
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FIG. 7: We quantify here the number of standard deviations 
n a for which the value of the form factor of the KIC deviates 
from the RMT prediction. This calculations were performed 
for I, 2, 3 and 4 kicks. The prominent feature is that for 
one kick we are always near 2 standard deviations away from 
the predicted result, systematically always undershooting the 
RMT result. For higher number of kicks the behavior is sta- 
tistically as expected. 

V(^2(l/r H ) 2 }coE = {K 2 (l/ m )y cOE = 2/N+0{N-*). 
Actually, in the limit Af = 00, few higher moments can 
be computed as well and COE distribution of ^(I/th) 
is conjectured to be exponential. In fig. [7] we plot the 
relative deviation of K^/th), for t = 1,2, 3, 4, from ex- 



pected COE average in terms of the number of expected 
standard deviations. It is clear that, at least for one 
kick, the deviation is exceeding COE model significantly. 
Namely we find the deviation in the same direction for 
all different numbers of qubits L, and for almost all L it 
is exceeding two standard deviations. We also find statis- 
tically significant deviations from COE for other values 
of t, in particular for t = 3, while the deviations for even 
arguments t = 2, 4 are less clear and conclusive. 

A specialized reader may inquire for a comparison with 
the DODO random matrix ensemble (which resembles 
semi-separable systems 10]). We note that the devia- 
tions from the COE observed in this article cannot be 
accounted by the semi-separable structure of the Floquet 
operator ((3|). 



V. CONCLUSIONS 

We have performed numerical calculations of large 
quasi-energy spectra of an interacting multi-qubit sys- 
tem, namely the kicked Ising chain. No analytical solu- 
tion of the model is known, i.e. the model is believed to 
be non-integrable. Consistently with previous results in 
the literature [ll| , we find good agreement of short-range 
level statistics of the model with Dyson's ensemble of cir- 
cular random matrices. However, when looking in detail 
at certain long-range spectral statistics, corresponding 
to short-times, we find notable and significant deviations 
from random matrix theory. This result reminds of non- 
universal regimes in semi-classical chaos widely studied 
in the 1990's. However, this behaviour cannot be at- 
tributed to periodic-orbits, since the system lacks any 
sensible definition of a classical limit. 

We believe that the numerical results are intriguing 
and await for theoretical explanation, perhaps in the di- 
rection of suggesting a new, abstract semi-classical pic- 
ture (perhaps along the lines of Ref.[r3|). 
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APPENDIX A: DIMENSIONS OF THE 
INVARIANT SUBSPACES 



Consider the computational basis S = 
{\rriomi . . . tul-x), rrij £ {0, 1, . . . , d — 1}} of the Hilbert 
space of L qudits H = H®^ dit . Let us generalize the 
translation operator allowing the rrij's to have integer 
values between and d — 1. The Hilbert space H is 
foliated into L subspaces Hk such that for any \ip) eHi, 

T\ip) = exp(27rifc/£)|V>) and H = 0fc =1 H fc . Let 
Pk be the orthogonal projection operator, such that 
Pk\ip) € Hk, for any H') € H. An elegant solution to 
the problem of calculating dim Hk is presented here, 
following (l3| . 

We first study the condition under which a state \n) £ 
S is projected to the null ket (zero vector). Let J be the 
smallest positive integer such that T J \n) = \n); we call J 
the primitive period of \n). The action of the projection 
operator Pk on \n) is 

P k \n) = ( J2^J,k) j ) \J2^,kT j ) |n), (Al) 



with (pi^k := e T . Notice that Pk\n) = if and only if 
4—1 

7 = 2JjL {tpj.k) 3 = 0. Since 7 is the sum of a geometric 
series, its calculation is straightforward: 7 = if and only 
if tpj,k 7^ 1. As a conclusion we obtain that Pk\n) 7^ if 
and only if ^ £ Z. 

Define the equivalence relation ~ in S as: \n) ~ \m) if 
there exists an integer j such that |n) = T J \m). Further- 
more, if |n) ~ |m), then Pu\n) oc Pj.|m), but if |n) ^ |m), 
then ((n\Pl)(Pk\m)) = 0. In other words, elements in 
different equivalence classes, are projected to orthogonal 



states. Thus counting the equivalence classes which are 
not projected to yields dimHk- 

Let N(J) be the number of equivalence classes whose 
elements have given primitive period J. If we call N(J) 
the number of elements in S that have primitive period 
J, then, N(J) = JN(J). Notice that 

£ N(J)= Yl JN(J)=d L (A2) 
{j|i£N} {Jl^eN} 

as the only allowed values for J are the divisors of L. 
Using Mobius inversion formula we obtain 

AT(J) = i £ M (£)<f" (A3) 

{m|ieN} 

Mobius function \x is defined over the positive integers as 
/i(l) = 1, = if n is divisible by the square of a 

prime, and in any other case, (J,(n) — (— l) p where p is 
the number of prime factors of n. Then, collecting our 
results, 

dimW fe = ^( J ) ( A4 ) 

since the only possible primitive periods J are the divi- 
sors of L. However the value of dim 7^. is well approxi- 
mated by 2 L /L for large values of L. 

APPENDIX B: THE OPTIMAL BASIS FOR 
DIAGONALIZATION OF KIC 

To get the spectra used in this paper, it is crucial to 
develop an optimal diagonalization scheme. Though the 
techniques relaying on Lanczos method [3] are fast they 
are not completely reliable. They loose precision as soon 
as some eigenvalues are close enough. In our experience 
the Lanczos method allows to obtain the full spectra for 
systems of up to 21 qubits, but the intrinsic numerical er- 
ror becomes comparable to the mean level spacing. Even 
for 18 qubits, the biggest numerical error in one level is 
already bigger than the smallest inter-level spacing. As 
we are performing very precise tests we require that our 
levels are highly reliable, making Lanczos a prohibitively 
inexact method. We prefer using direct diagonalization 
with specialized routines [l5j . 

Let U StK i = U kick {b/2)Ui s i ng (J)U kick (b/2) be the sym- 
metrized version of Uki, which however has the same 
spectrum due to unitary equivalence. Using an appro- 
priate basis is important both to take advantage of the 
natural block diagonal decomposition of t/ Sj Ki (due to 
symmetry P) and its symmetric character (due to sym- 
metry JC). The basis is constructed as follows. Let \n) 
be a representative of a class defined by ~ (see App. [A"f. 
and such that P&|n) ^ 0. If (n\RPk\n) = then both 
Pk\n) ± JCRPk\n) are used as members of the basis. If 
(n\RPk\n) at least one of Pk\n) ± JCRPk\n) is not 
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the null kct, and can be incorporated into the base. By plete basis in Hk- Moreover, this basis is invariant under 
choosing \n) from all different classes we build a com- K, and is orthogonal. 



